## ============================================================================
## 06_si_s2_audience_effects.R
## ----------------------------------------------------------------------------
## Purpose:   Supplementary analyses for Study 1 (audience effects):
##            formal ATE tables, equivalence tests, covariate-adjusted
##            estimates, pooled estimator, and ATE-to-partisan-gap
##            benchmarking. Corresponds to SI Section S2.2.
##
## Requires:  01_manuscript.R must be run first (reads ATEs CSV).
##
## Inputs:    combined_data.rds (via 00_setup.R),
##            ates_study1_by_sample_unadj.rds (from 01_manuscript.R)
## Outputs:   Tables/table_s23.tex,
##            Tables/table_s24.tex,
##            Tables/table_s25.tex,
##            Tables/table_s26.tex,
##            Figures/fig_s12.pdf,
##            ates_study1_pooled_unadj.rds
## ============================================================================

source("00_setup.R")

## ---- Table S23: Audience effects by sample + between-county differences ----
## Reconstruct unadjusted ATEs from CSV (written by 01_manuscript.R) 
summary_df <- read_rds("ates_study1_by_sample_unadj.rds")

## Adjust CIs and P-Values for multiple comparisons
summary_df_unadj <-
  summary_df %>%
  filter(outcome_group != "Indigenous Awareness") %>%
  mutate(
    ## Adjust CIs
    ci.lwr95_adj = case_when(
      sample == "Australia" ~ estimate - qnorm(1 - (0.05 / 7) / 2) * std.error,
      sample == "United States" ~ estimate - qnorm(1 - (0.05 / 6) / 2) * std.error,
    ),
    ci.upr95_adj = case_when(
      sample == "Australia" ~ estimate + qnorm(1 - (0.05 / 7) / 2) * std.error,
      sample == "United States" ~ estimate + qnorm(1 - (0.05 / 6) / 2) * std.error,
    ),
    ## Adjust CIs
    ci.lwr90_adj = case_when(
      sample == "Australia" ~ estimate - qnorm(1 - (0.10 / 7) / 2) * std.error,
      sample == "United States" ~ estimate - qnorm(1 - (0.10 / 6) / 2) * std.error,
    ),
    ci.upr90_adj = case_when(
      sample == "Australia" ~ estimate + qnorm(1 - (0.10 / 7) / 2) * std.error,
      sample == "United States" ~ estimate + qnorm(1 - (0.10 / 6) / 2) * std.error,
    )
  ) %>%
  group_by(sample) %>%
  mutate(## Adjust P-Values
    pval_adj = case_when(
      sample == "Australia" ~ p.adjust(p.value, method = "bonferroni"),
      sample == "United States" ~ p.adjust(p.value, method = "bonferroni")
    )) %>%
  ungroup() %>%
  select(
    outcome_group,
    sample,
    estimate,
    std.error,
    p.value,
    ci.lwr95_adj,
    ci.upr95_adj,
    pval_adj,
    ci.lwr90_adj,
    ci.upr90_adj
  ) %>%
  arrange(outcome_group)

## Calculated between country differences for SI
diffs_df <-
  summary_df_unadj %>%
  select(outcome_group, sample, estimate, std.error) %>%
  filter(outcome_group != "Voice Referendum") %>%
  pivot_wider(
    names_from = sample,
    values_from = c(estimate, std.error),
    names_sep = "_"
  ) %>%
  mutate(
    est_diff = `estimate_Australia` - `estimate_United States`,
    se_diff = sqrt(`std.error_Australia`^2 + `std.error_United States`^2),
    p_diff = 2 * pnorm(-abs(est_diff / se_diff)),
    ci_lwr_diff = est_diff - qnorm(0.975) * se_diff,
    ci_upr_diff = est_diff + qnorm(0.975) * se_diff,

    ## Bonferroni simultaneous 95% CIs (FWER control)
    ci_lwr_diff_adj = est_diff - qnorm(1 - (0.05 / 6) / 2) * se_diff,
    ci_upr_diff_adj = est_diff + qnorm(1 - (0.05 / 6) / 2) * se_diff,
    p_diff_adj = p.adjust(p_diff, method = 'bonferroni'),
    sample = "Difference"
  ) %>%
  select(outcome_group, contains("_diff"))

est_diff_tab <-
  summary_df_unadj %>%
  mutate(entry = table_entry(est = estimate, se = std.error)) %>%
  select(outcome_group, sample, entry) %>%
  pivot_wider(names_from = sample, values_from = entry, values_fill = "-") %>%
  left_join(
    diffs_df %>%
      mutate(
        diff = table_entry(est = est_diff, se = se_diff),
        diff_ci = paste0("[", round(ci_lwr_diff_adj, 2), ", ", round(ci_upr_diff_adj, 2), "]"),
        diff_p = case_when(
          p_diff_adj < 0.001 ~ "< 0.001",
          round(p_diff_adj, 3) == 0.001 ~ "0.001",
          p_diff_adj > 0.99 ~ "> 0.99",
          .default = sprintf("%.3f", p_diff_adj)
        )
      ) %>%
      select(outcome_group, diff, diff_ci, diff_p),
    by = "outcome_group"
  ) %>%
  arrange(desc(outcome_group))

## Export for latex
est_diff_tab %>%
  mutate(
    outcome_group = paste0("\\textit{", outcome_group, "}"),
    across(where(is.character), ~replace_na(., "-"))  # Voice Referendum: no US/diff cols
  ) %>%
  xtable() %>%
  print(
    include.rownames = FALSE,
    include.colnames = FALSE,
    only.contents = TRUE,
    type = "latex",
    sanitize.text.function = identity,
    comment = FALSE,
    hline.after = NULL,
    file = paste0(tab_dir, "table_s23.tex")
  )

## ---- Table S24: Benchmarking ATEs against partisan gaps -----

## Subset to untreated respondents to calculate partisan gaps
subset_dat <-
  combined_dat %>%
  filter(Z_land_info == 0)  %>%
  mutate(sample = factor(
    sample,
    levels = c("AU", "US"),
    labels = c("Australia", "United States")
  ))

pid_s_vars <- c("y_substantive_idx_s", "y_symbolic_idx_s", "y_guilt_idx_s",
                "y_cbr_idx_s", "y_acknowledge_s", "y_infoseek_s",
                "y_recall_binary_s", "y_voice_vote_s")

dim_pid_s <- map_dfr(pid_s_vars, function(y_var) {
  dat <- if (y_var == "y_voice_vote_s") {
    subset_dat %>% filter(sample == "Australia")
  } else {
    subset_dat
  }
  dat %>%
    group_by(x_pid_comb2, sample) %>%
    group_modify(~ tidy(lm_robust(reformulate("1", y_var), data = .x))) %>%
    ungroup() %>%
    filter(term == "(Intercept)") %>%
    mutate(outcome = y_var)
}) %>%
  mutate(outcome = factor(
    outcome,
    levels = c("y_recall_binary_s", "y_acknowledge_s", "y_infoseek_s",
               "y_guilt_idx_s", "y_cbr_idx_s", "y_symbolic_idx_s",
               "y_substantive_idx_s", "y_voice_vote_s"),
    labels = c(
      "Indigenous Awareness",
      "Indigenous Land",
      "Information Seeking",
      "Collective Guilt",
      "Colorblind Racism",
      "Symbolic Reparations",
      "Substantive Reparations",
      "Voice Referendum"
    )
  ),
  x_pid = factor(x_pid_comb2, levels = c("Left", "Independent", "Right"))
  )


pid_diff_s <-
  dim_pid_s %>%
  filter(x_pid != "Independent") %>%
  select(x_pid, sample, outcome, estimate, std.error, p.value) %>%
  pivot_wider(names_from = x_pid, values_from = c(estimate, std.error, p.value)) %>%
  mutate(
    # Calculate between country differences for each outcome
    diff_estimate = `estimate_Left` - `estimate_Right`,
    diff_std_error = sqrt(`std.error_Left`^2 + `std.error_Right`^2),
    diff_p_value = 2 * pnorm(abs(diff_estimate) / diff_std_error, lower.tail = FALSE) # Two-tailed p-value
  )

pid_diff_s %>%
  select(outcome, sample, starts_with("diff_")) %>%
  arrange(outcome)

## Combined with ATEs
ate_dat <- read_rds("ates_study1_by_sample_unadj.rds")

benchmark_dat <-
  ate_dat %>%
  select(outcome_group, sample, estimate) %>%
  rename(outcome = outcome_group) %>%
  left_join(
    pid_diff_s %>%
      select(outcome, sample, diff_estimate),
    by = c("outcome", "sample")
  ) %>%
  rename(ate = estimate,
         pid_diff = diff_estimate) %>%
  mutate(
    ratio = abs(ate/pid_diff),
    outcome = factor(
      outcome,
      levels = c(
        "Indigenous Awareness",
        "Indigenous Land",
        "Information Seeking",
        "Collective Guilt",
        "Colorblind Racism",
        "Symbolic Reparations",
        "Substantive Reparations",
        "Voice Referendum"
      )
    ))

latex_tab_pid <-
  benchmark_dat %>%
  filter(!(outcome %in% c("Indigenous Awareness"))) %>%
  arrange(sample, outcome) %>%
  mutate(across(where(is.numeric), \(x) ifelse(is.na(x), NA_character_,
                                               sprintf("%.2f", round(x, 2))))) %>%
  pivot_wider(names_from = sample, values_from = c(ate, pid_diff, ratio),
              values_fill = "-") %>%
  select(outcome, ends_with("_Australia"), ends_with("_United States"))


latex_tab_pid %>%
  mutate(
    outcome = paste0("\\textit{", outcome, "}"),
    across(where(is.character), ~replace_na(., "-"))  # Voice Referendum: no US cols
  ) %>%
  xtable() %>%
  print(
    include.rownames = FALSE,
    include.colnames = FALSE,
    only.contents = TRUE,
    type = "latex",
    sanitize.text.function = identity,
    comment = FALSE,
    hline.after = NULL,
    file = paste0(tab_dir, "table_s24.tex")
  )

## ---- Estimate covariate-adjusted ATEs via crossEstimation for Table S25 ----

## NB: geographic variables are x_state_comb for Australian sample and
##     x_region for us sample
us_covs <-
  c(
    "x_edu_college",
    "x_employ_binary",
    "x_gender_f",
    "x_age_cat",
    "x_hhi",
    "x_native",
    "x_pid_comb2",
    "x_region",
    "x_ethnicity"
  )

au_covs <-
  c(
    "x_edu_college",
    "x_employ_binary",
    "x_gender_f",
    "x_age_cat",
    "x_hhi",
    "x_native",
    "x_pid_comb2",
    "x_state_comb",
    "x_ideo_cat"
  )

au_outcomes <- c(
  "y_guilt_idx_s",
  "y_cbr_idx_s",
  "y_symbolic_idx_s",
  "y_substantive_idx_s",
  "y_voice_vote_s",
  "y_infoseek_s",
  "y_acknowledge_s"
)

us_outcomes <- c(
  "y_guilt_idx_s",
  "y_cbr_idx_s",
  "y_symbolic_idx_s",
  "y_substantive_idx_s",
  "y_infoseek_s",
  "y_acknowledge_s"
)

est_us_df <-
  combined_dat %>%
  filter(sample == "US") %>%
  select(
    all_of(us_outcomes),
    all_of(us_covs),
    Z_land_info
  ) %>%
  filter(!is.na(Z_land_info))

est_au_df <-
  combined_dat %>%
  filter(sample == "AU") %>%
  select(
    all_of(au_outcomes),
    all_of(au_covs),
    Z_land_info
  ) %>%
  filter(!is.na(Z_land_info))


## Loop over all the outcome measures, apply the algorithm and store the
## results.
set.seed(123)

au_forest_est <- list()
us_forest_est <- list()

for (i in 1:length(au_outcomes)) {

  ## Filter out all complete cases for each outcome
  est_au_adj <-
    est_au_df %>%
    filter(if_all(c(au_outcomes[i], all_of(au_covs)), \(x) complete.cases(x)))

  X <-
    model.matrix(as.formula(paste("~", paste(
      au_covs, collapse = "+"
    ))),
    data = est_au_adj)

  ## Specify outcome for this loop
  y <- est_au_adj %>% pull(au_outcomes[i]) %>% as.numeric()

  ## Binary treatment indicator
  W <- as.numeric(est_au_adj$Z_land_info)

  tmp <-
    suppressWarnings(ate.randomForest(
      X = X,
      Y = y,
      W = W,
      conf.level = 0.95
    ))

  au_forest_est[[i]] <-
    data.frame(
      estimate = tmp$tau,
      std.error = sqrt(tmp$var),
      conf.low = tmp$conf.int[1],
      conf.high = tmp$conf.int[2],
      outcome = au_outcomes[i],
      sample = "Australia"
    )

}

for (i in 1:length(us_outcomes)) {

  ## Filter out all complete cases for each outcome
  est_us_adj <-
    est_us_df %>%
    filter(if_all(c(us_outcomes[i], all_of(us_covs)), \(x) complete.cases(x)))

  X <-
    model.matrix(as.formula(paste("~", paste(
      us_covs, collapse = "+"
    ))),
    data = est_us_adj)

  ## Specify outcome for this loop
  y <- est_us_adj %>% pull(us_outcomes[i]) %>% as.numeric()

  ## Binary treatment indicator
  W <- as.numeric(est_us_adj$Z_land_info)

  tmp <-
    suppressWarnings(ate.randomForest(
      X = X,
      Y = y,
      W = W,
      conf.level = 0.95
    ))

  us_forest_est[[i]] <-
    data.frame(
      estimate = tmp$tau,
      std.error = sqrt(tmp$var),
      conf.low = tmp$conf.int[1],
      conf.high = tmp$conf.int[2],
      outcome = us_outcomes[i],
      sample = "United States"
    )

}


## ---- Table S25: Estimates with and without covariate adjustment -----

## Combine estimates and export table of results
summary_df_adj <-
  bind_rows(au_forest_est) %>%
  bind_rows(us_forest_est) %>%
  rename(outcome_group = outcome) %>%
  mutate(
    estimator = "Adjusted",
    statistic = estimate / std.error,
    p.value = 2 * (1 - pnorm(abs(statistic))),
    outcome_group = case_when(
      outcome_group == "y_guilt_idx_s" ~ "Collective Guilt",
      outcome_group == "y_cbr_idx_s" ~ "Colorblind Racism",
      outcome_group == "y_symbolic_idx_s" ~ "Symbolic Reparations",
      outcome_group == "y_substantive_idx_s" ~ "Substantive Reparations",
      outcome_group == "y_voice_vote_s" ~ "Voting for Referendum",
      outcome_group == "y_infoseek_s" ~ "Information Seeking",
      outcome_group == "y_acknowledge_s" ~ "Indigenous Land"
    ),
    outcome_group = factor(
      outcome_group,
      levels = c(
        "Voting for Referendum",
        "Symbolic Reparations",
        "Substantive Reparations",
        "Colorblind Racism",
        "Collective Guilt",
        "Information Seeking",
        "Indigenous Land"
      ),
      labels = c(
        "Voice Referendum",
        "Symbolic Reparations",
        "Substantive Reparations",
        "Colorblind Racism",
        "Collective Guilt",
        "Information Seeking",
        "Indigenous Land"
      )
    ),
    sample = factor(
      sample,
      levels = c("Australia", "United States")
    ),
    ## Adjust CIs and P-Values for multiple comparisons
    ci.lwr95_adj = case_when(
      sample == "Australia" ~ estimate - qnorm(1 - (0.05 / 7) / 2) * std.error,
      sample == "United States" ~ estimate - qnorm(1 - (0.05 / 6) / 2) * std.error,
    ),
    ci.upr95_adj = case_when(
      sample == "Australia" ~ estimate + qnorm(1 - (0.05 / 7) / 2) * std.error,
      sample == "United States" ~ estimate + qnorm(1 - (0.05 / 6) / 2) * std.error,
    ),
    ## Adjust CIs
    ci.lwr90_adj = case_when(
      sample == "Australia" ~ estimate - qnorm(1 - (0.10 / 7) / 2) * std.error,
      sample == "United States" ~ estimate - qnorm(1 - (0.10 / 6) / 2) * std.error,
    ),
    ci.upr90_adj = case_when(
      sample == "Australia" ~ estimate + qnorm(1 - (0.10 / 7) / 2) * std.error,
      sample == "United States" ~ estimate + qnorm(1 - (0.10 / 6) / 2) * std.error,
    )
  ) %>%
  group_by(sample) %>% 
  mutate(## Adjust P-Values
    pval_adj = case_when(
      sample == "Australia" ~ p.adjust(p.value, method = "bonferroni"),
      sample == "United States" ~ p.adjust(p.value, method = "bonferroni")
    )) %>%
  ungroup() %>%
  select(
    outcome_group,
    sample,
    estimator,
    estimate,
    std.error,
    p.value,
    ci.lwr95_adj,
    ci.upr95_adj,
    pval_adj,
    ci.lwr90_adj,
    ci.upr90_adj
  ) %>%
  arrange(outcome_group)


summary_combined <-
  summary_df_unadj %>%
  mutate(estimator = "Unadjusted") %>%
  bind_rows(summary_df_adj) 


### Export table of results
compare_tab <- 
  summary_combined %>%
  mutate(
    entry = table_entry(est = estimate, se = std.error),
    pval = case_when(
      pval_adj < 0.001 ~ "< 0.001",
      round(pval_adj, 3) == 0.001 ~ "0.001",
      pval_adj > 0.99 ~ "> 0.99",
      .default = sprintf("%.3f", pval_adj)
    )
  ) %>%
  select(sample, outcome_group, estimator, entry, pval) %>%
  pivot_wider(names_from = c(estimator, sample), values_from = c(entry, pval),
              values_fill = "-") %>%
  arrange(desc(outcome_group)) %>%
  mutate(
    outcome_group = paste0("\\textit{", outcome_group, "}"),
    across(where(is.character), ~replace_na(., "-"))  # Voice Referendum: no US cols
  ) 

compare_tab %>%
  select(
    outcome_group,
    entry_Unadjusted_Australia,
    entry_Adjusted_Australia,   
    `entry_Unadjusted_United States`, 
    `entry_Adjusted_United States`
  ) %>% 
  xtable() %>%
  print(
    include.rownames = FALSE,
    include.colnames = FALSE,
    only.contents = TRUE,
    type = "latex",
    sanitize.text.function = identity,
    comment = FALSE,
    hline.after = NULL,
    file = paste0(tab_dir, "table_s25.tex")
  )


## ---- Table S26: Audience effects using pooled estimator -----

## Re-standardize outcomes across the combined (pooled) sample
est_df_pooled <-
  combined_dat %>%
  mutate(
    y_guilt_idx_s = glass_delta(y_guilt_idx, Z = Z_land_info, reference = 0),
    y_cbr_idx_s = glass_delta(y_cbr_idx, Z = Z_land_info, reference = 0),
    y_substantive_idx_s = glass_delta(y_substantive_idx, Z = Z_land_info, reference = 0),
    y_symbolic_idx_s = glass_delta(y_symbolic_idx, Z = Z_land_info, reference = 0),
    y_voice_vote_s = glass_delta(y_voice_vote, Z = Z_land_info, reference = 0),
    y_infoseek_s = glass_delta(y_infoseek_n, Z = Z_land_info, reference = 0),
    y_acknowledge_s = glass_delta(y_acknowledge_n, Z = Z_land_info, reference = 0),
    y_recall_binary_s = glass_delta(y_recall_binary, Z = Z_land_info, reference = 0),
  ) %>%
  select(
    sample,
    y_guilt_idx_s,
    y_cbr_idx_s,
    y_substantive_idx_s,
    y_symbolic_idx_s,
    y_voice_vote_s,
    y_infoseek_s,
    y_recall_binary_s,
    y_acknowledge_s,
    Z_land_info
  ) %>%
  pivot_longer(cols = starts_with("y_"),
               names_to = "outcome_group",
               values_to = "y") %>%
  mutate(
    Z = ifelse(Z_land_info == 1, "Treatment", "Control"),
    outcome_group = case_when(
      outcome_group == "y_guilt_idx_s" ~ "Collective Guilt",
      outcome_group == "y_cbr_idx_s" ~ "Colorblind Racism",
      outcome_group == "y_symbolic_idx_s" ~ "Symbolic Reparations",
      outcome_group == "y_substantive_idx_s" ~ "Substantive Reparations",
      outcome_group == "y_voice_vote_s" ~ "Voting for Referendum",
      outcome_group == "y_infoseek_s" ~ "Information Seeking",
      outcome_group == "y_acknowledge_s" ~ "Indigenous Land",
      outcome_group == "y_recall_binary_s" ~ "Indigenous Awareness"
    )
  ) %>%
  filter(!is.na(Z))

tmp_pooled <-
  est_df_pooled %>%
  filter(outcome_group == "Voting for Referendum", sample == "AU") %>%
  {tidy(lm_robust(y ~ Z, data = .))} %>%
  filter(term != "(Intercept)") %>%
  mutate(outcome_group = "Voting for Referendum")

summary_df_pooled_raw <-
  est_df_pooled %>%
  filter(outcome_group != "Voting for Referendum") %>%
  group_by(outcome_group) %>%
  group_modify(~ tidy(lm_lin(
    y ~ Z, covariates = ~ sample, data = .x
  ))) %>%
  filter(term == "ZTreatment") %>%
  bind_rows(tmp_pooled) %>%
  mutate(outcome_group = factor(
    outcome_group,
    levels = c(
      "Voting for Referendum",
      "Substantive Reparations",
      "Symbolic Reparations",
      "Colorblind Racism",
      "Collective Guilt",
      "Information Seeking",
      "Indigenous Land",
      "Indigenous Awareness"
    ),
    labels = c(
      "Voice Referendum",
      "Substantive Reparations",
      "Symbolic Reparations",
      "Colorblind Racism",
      "Collective Guilt",
      "Information Seeking",
      "Indigenous Land",
      "Indigenous Awareness"
    )
  )) %>%
  ungroup()

## Adjust CIs and P-Values for multiple comparisons
summary_df_pooled <-
  summary_df_pooled_raw %>%
  mutate(
    ## Adjust CIs
    ci.lwr95_adj = estimate - qnorm(1 - (0.05 / 7) / 2) * std.error,
    ci.upr95_adj = estimate + qnorm(1 - (0.05 / 7) / 2) * std.error,
    ## Adjust CIs
    ci.lwr90_adj = estimate - qnorm(1 - (0.10 / 7) / 2) * std.error,
    ci.upr90_adj = estimate + qnorm(1 - (0.10 / 7) / 2) * std.error,
    pval_adj = p.adjust(p.value, method = "bonferroni")
  ) %>%
  mutate(across(where(is.numeric), \(x) round(x, 3))) %>%
  select(
    outcome_group,
    estimate,
    std.error,
    p.value,
    ci.lwr95_adj,
    ci.upr95_adj,
    pval_adj,
    ci.lwr90_adj,
    ci.upr90_adj
  ) %>%
  arrange(outcome_group)

write_rds(summary_df_pooled, "ates_study1_pooled_unadj.rds")

latex_tab_pooled <-
  summary_df_pooled %>%
  filter(outcome_group != "Indigenous Awareness") %>%
  arrange(desc(outcome_group)) %>%
  mutate(
    outcome_group = paste0("\\textit{", outcome_group, "}"),  # Italicize text
    entry = table_entry(est = estimate, se = std.error),
    ci95 = paste0("[", round(ci.lwr95_adj, 2), ", ", round(ci.upr95_adj, 2), "]"),
    ci90 = paste0("[", round(ci.lwr90_adj, 2), ", ", round(ci.upr90_adj, 2), "]"),
    pval = case_when(
      pval_adj < 0.001 ~ "< 0.001",
      round(pval_adj, 3) == 0.001 ~ "0.001",
      pval_adj > 0.99 ~ "> 0.99",
      .default = sprintf("%.3f", pval_adj)
    )
  ) %>%
  select(outcome_group, entry, ci95, ci90, pval)

latex_tab_pooled %>%
  xtable() %>%
  print(
    include.rownames = FALSE,
    include.colnames = FALSE,
    only.contents = TRUE,
    type = "latex",
    sanitize.text.function = identity,
    comment = FALSE,
    hline.after = NULL,
    file = paste0(tab_dir, "table_s26.tex")
  )

## ---- Figure S12: Audience effects using pooled estimator ------
p_pooled <-
  summary_df_pooled %>%
  filter(outcome_group != "Indigenous Awareness") %>%
  ggplot(data = ., aes(x = estimate, y = outcome_group)) +
  geom_vline(
    xintercept = 0,
    color = "black",
    linetype = "solid",
    linewidth = 0.2
  ) +
  geom_vline(
    xintercept = -0.2,
    color = "black",
    linetype = "dotted",
    linewidth = 0.2
  ) +
  geom_vline(
    xintercept = 0.2,
    color = "black",
    linetype = "dotted",
    linewidth = 0.2
  ) +
  geom_errorbar(
    aes(xmin = ci.lwr95_adj, xmax = ci.upr95_adj),
    width = 0,
    linewidth = .75,
    orientation = "y"
  ) +
  geom_errorbar(
    aes(xmin = ci.lwr90_adj, xmax = ci.upr90_adj),
    width = 0,
    linewidth = 1.5,
    orientation = "y"
  ) +
  geom_point(
    size = 4,
    fill = "black",
    pch = 22,
    color = "white"
  ) +
  scale_x_continuous(limits = c(-0.32, 0.32),
                     breaks = c(-0.30, -0.20, -0.10, 0, 0.10, 0.20, 0.30)) +
  theme_tufte(base_size = 14) +
  labs(subtitle = "",
       x = "Average treatment effect (in standard units)",
       y = "") +
  theme(
    axis.ticks.y = element_blank(),
    axis.ticks.x = element_line(linewidth = 0.25),
    axis.line.x = element_line(linewidth = 0.25),
    axis.text.x = element_text(size = 12, color = "black"),
    axis.text.y = element_text(size = 12, color = "black"),
    axis.title.x = element_text(size = 12, color = "black"),
    strip.text = element_text(size = 14, hjust = 0.5),
    plot.margin = unit(c(
      t = -1,
      r = 0.2,
      b = 0,
      l = -1
    ), "lines"),
    legend.position = "none",
    legend.margin = margin(0, 0, 0, 0),
    legend.box.spacing = unit(0, "pt"),
    legend.box.margin = margin(0, 0, 0, 0)
  )

p_pooled

ggsave(
  filename = paste0(fig_dir, "fig_s12.pdf"),
  p_pooled,
  width = 5.5,
  height = 4.5
)

## ---- Session info -----
sessionInfo()
